function f = func_Rb_int(rbv_sol,kv_st,nv_st,cv_st,hv_st,a_st,...
                         beta,lambda,gamma,delta,rho_a,sigma_a,alpha,psi,nu,z_trapz,n_trapz,expmx,eta1_c,eta2_c,eta3_c)


        % numerical integration to compute Rbar explicitly
        rhsee_each = zeros(n_trapz,1);
        for z=1:n_trapz
        
            a_next  = rho_a*a_st + sigma_a*z_trapz(z);
            h_next  = ((1-alpha)/psi*exp(a_next)*kv_st^(alpha))^(1/(alpha+1/nu));
            y_next  = exp(a_next) * kv_st^(alpha)*h_next^(1-alpha);
        
            rkhats_next = log ( rbv_sol * (1-1/(kv_st/nv_st)) * ( 1+lambda*(1-gamma) ) - (1-delta) );
            rkhat_next  = log ( alpha*((1-alpha)/psi)^((1-alpha)/(alpha+1/nu)) ) ...
                        + (1+nu)/(alpha*nu+1)*a_next - (1-alpha)/(alpha*nu+1)*log(kv_st);
            xi_next     = exp(rkhat_next)/exp(rkhats_next);
            
            K  = [kv_st rbv_sol nv_st xi_next];
            x0 = (log(K)).^expmx;
            x2_next  = prod(x0,2)';
            
            if xi_next >= 1
                                
                pibv_next  = ( exp(rkhat_next) + 1-delta )*(kv_st/nv_st)*nv_st - rbv_sol*((kv_st/nv_st)-1)*nv_st - nv_st;
                
                if pibv_next > 0
                    c_next     = exp(x2_next*eta1_c);
%    !                rv_next    = rk_next * lv(t)/(lv(t)-1) - rbv(t)*lambda
                    rhsee_temp = beta / ( c_next - psi*h_next^(1+1/nu)/(1+1/nu) ) * rbv_sol;
                else
                    c_next     = exp(x2_next*eta3_c);
%    !                rv_next    = rk_next * lv(t)/(lv(t)-1) - rbv(t)*lambda
                    rhsee_temp = beta / ( c_next - psi*h_next^(1+1/nu)/(1+1/nu) ) * rbv_sol;
                end    

            else % if run happens
                
                c_next     = exp(x2_next*eta2_c);
                rv_next    = ( exp(rkhat_next) + 1-delta ) * (kv_st/nv_st)/((kv_st/nv_st)-1) - rbv_sol*lambda;
                rhsee_temp = beta / ( c_next - psi*h_next^(1+1/nu)/(1+1/nu) ) * rv_next;
                
            end

            % multiply density, and cut off the both tails outside of z_trapz
            pdf = normpdf(z_trapz(z));
            cdf = normcdf(z_trapz(1));
            rhsee_each(z) = rhsee_temp * pdf / ( 1 - 2*cdf );

        end
        
        % compute integeral
        rhsee_derived = trapz(z_trapz,rhsee_each);
        f =  1/( cv_st - psi*hv_st^(1+1/nu)/(1+1/nu) ) - rhsee_derived;

end